double  interaccion(caja *bx,double rc)
{
	int ind;
	int ank, ank2, ankh;
	ank=6;
	ank2=2*ank;
	ankh=ank/2;
	double dist2, dist,rqInv,dr[3], tsn, ccel3,potdm, ePoten;
	double conrep= 4.0*pow(1.0,ank2); //Parte repulsiva del potencial Leonard-Jones
	double conref= ank2*conrep; //Parte repulsiva de la fuerza del Leonard-Jones
	double conat= 4.0*pow(1.0,ank); //Parte atractiva del potencial Leonard-Jones
	double conatf= ank*conat;// Parte atractiva de la fuerza del Leonard-Jones

	bx->initForces();

	double LJcut= pow(1.0/(rc*rc),ankh)*(pow(1.0/(rc*rc),ankh)*conrep-conat);

	ePoten=0.;
	for(int i = 0;i<bx->npart;i++)
	//for(int i = 0;i<bx->npart-1;i++)
	{
		//cout<<i<<"  "<<bx->part[i].nvecinos<<endl;
		for(int j=0;j<bx->part[i].nvecinos;j++)
		//for(int j=i+1;j<bx->npart;j++)
		{

			ind = bx->part[i].lstVecinos[j];
			//cout<<i<<"  **** "<<ind<<endl;
			//ind = j;
			dist2= bx->distancia2(i,ind);
			//if(dist2<0.54) dist2=0.54;
			rqInv =1.0/dist2;//1/dist^2
			tsn=pow(rqInv,ankh);
			//F/r=-(dPOTDM/dr)/r (CORRECTO)
			ccel3 =rqInv*tsn*(conref*tsn-conatf);
			//ccel3=0;
			//cout<<dist2<<endl;
			potdm = tsn*(conrep*tsn-conat);
//			if(i==2 )
//			{
////				cout<<"conrep: "<<conrep<<endl;
////				cout<<"conat: "<<conat<<endl;
////				cout<<"dis2: "<<dist2<<endl;
//				cout<<"****** "<<potdm<<"  "<<tsn<<"  "<<sqrt(dist2)<<endl;
//				//system("PAUSE");
//			}


			ePoten+=potdm;
			for (int k=0;k<3;k++)
			{
				dr[k]= bx->dr[k];
				//Para la particula I
				bx->part[i].force[k]+=ccel3*dr[k];
				//Para la particula J
				bx->part[j].force[k]-=ccel3*dr[k];
			}
			//energia(dist);
		}
	}
	ePoten-=LJcut*(bx->npart)*((bx->npart)-1)/2.0;
	return ePoten;

}


